Unique Transcriptional Signatures Correlate with Behavioral and Psychological Symptom Domains in Alzheimer’s Disease

Despite the significant burden, cost, and worse prognosis of Alzheimer’s disease (AD) with behavioral and psychological symptoms of dementia (BPSD), little is known about the molecular causes of these symptoms. Using antemortem assessments of BPSD in AD, we demonstrate that individual BPSD can be grouped into 4 domain factors in our sample: affective, apathy, agitation, and psychosis. Then, we performed a transcriptome-wide analysis for each domain utilizing bulk RNA-seq of post-mortem anterior cingulate cortex (ACC) tissue. Though all 4 domains are associated with a predominantly downregulated pattern of hundreds of differentially expressed genes (DEGs), most DEGs are unique to each domain, with only 22 DEGs being common to all BPSD domains, including TIMP1. Weighted gene co-expression network analysis (WGCNA) yielded multiple transcriptional modules that were shared between BPSD domains or unique to each domain, and NetDecoder was used to analyze context-dependent information flow through the biological network. For the agitation domain, we found that all DEGs and a highly correlated transcriptional module were functionally enriched for ECM-related genes including TIMP1, TAGLN, and FLNA. Another unique transcriptional module also associated with the agitation domain was enriched with genes involved in post-synaptic signaling, including DRD1, PDE1B, CAMK4, and GABRA4. By comparing context-dependent changes in DEGs between cases and control networks, ESR1 and PARK2 were implicated as two high impact genes associated with agitation that mediated significant information flow through the biological network. Overall, our work establishes unique targets for future study of the biological mechanisms of BPSD and resultant drug development.


INTRODUCTION
Although dementia is defined largely by memory and cognitive decline which results in the loss of the ability to function independently, behavioral and psychological symptoms in dementia (BPSD) play a large role in a patient's overall functional level [1]. Significant BPSD are more the rule than the exception, as >90% of people with dementia develop BPSD during the disease course [2,3]. BPSD encompass a wide array of symptoms and include aggression, agitation, hyperactivity, compulsions, disinhibition, anxiety, depression and dysphoria, euphoria delusions, and hallucinations. In addition, there is evidence that mild behavioral impairment (MBI) is the first sign of dementia in some people, analogous to the more widely recognized mild cognitive impairment (MCI) [4,5]. In agreement with this, neuropsychiatric symptoms can precede a dementia diagnosis [4,6], with some estimates suggesting that over half demonstrate neuropsychiatric symptoms before a diagnosis of a cognitive disorder, including MCI [7]. Neuropsychiatric symptoms are also associated with faster cognitive decline in cognitively unimpaired individuals [8][9][10] and those with MCI [11], quicker development of dementia in people with MCI [12][13][14], and faster dementia progression [15]. In addition, BPSD result in poorer quality of life for those with dementia [16,17] as well as significant caregiver distress, often greater than with cognitive deficits [18,19]. Emergent and difficult to treat, BPSD are often the cause of hospitalization and institutionalization for persons with dementia [20,21]. Behavioral interventions remain first-line for treating any BPSD, and there is still no FDA approved medications that are indicated to treat any BPSD [22].
BPSD are a heterogeneous group of symptoms, and while each symptom has its own frequency of presentation, studies have suggested that certain symptoms co-occur at greater rates than others [23,24]. In particular, a systematic review of 62 studies utilizing unbiased clustering of BPSD generally found that affective symptoms (dysphoria, anxiety), apathy, hyperactivity-impulsivity-disinhibition-agitation-aggression (HIDA), and psychosis (delusions, hallucinations) all tend to form independent clusters [24]. This could implicate similar molecular mechanisms underlying each cluster. Despite the ubiquity and burden of BPSD clinically, very few molecular studies exist to elucidate the molecular mechanisms underlying BPSD [25][26][27], with the exception of AD with psychosis [28][29][30]. To begin to identify the molecular mechanisms of BPSD, in this study, we first verified the occurrence of BPSD domains in a sample of older adults based on the use of a structured clinical interview administered within two years of death. Then, we confirmed clustering of BPSD into four predominant domains and then performed bulk RNAseq analysis on the post-mortem anterior cingulate cortex (ACC) from a subset of those with AD and varying burdens of BPSD. Finally, we performed weighted gene co-expression network analysis (WGCNA) to identify transcriptional clusters associated with each BPSD domain. Finally, we used an algorithm to assess contextdependent information flow differences between the transcriptional networks of cases and controls to determine likely molecular drivers behind each BPSD domain. Our data set to better understand these four BPSD domains on the transcriptional level and yield promising targets for future mechanistic studies and novel therapeutics.

Subjects
Community-dwelling older adults who later developed dementia and their informants were recruited through the Rush Alzheimer's Disease Center (RADC) memory clinic. As a part of their research visits, trained research assistants conducted standardized clinical interviews via telephone with informants, including assessment of the frequency and severity of numerous BPSD. Questions pertaining to BPSD were developed by two neuropsychologists from clinical descriptions in the literature, observations of patient behaviors, and caregiver interviews, as described previously [31][32][33][34].
This study was approved by an Institutional Review Board of Rush University Medical Center (RUMC). For all subjects, consent for brain autopsy was obtained after death from next of kin and a witness by RUMC staff.
Inclusion criteria for subjects in the present study were based on the NINCDS-ADRDA clinical criteria for "probable AD" with diagnosticians blinded to post-mortem findings [35]. Namely, all subjects had a history of cognitive decline, impairment in memory and at least one other cognitive domain, and no other conditions judged to be probably contributing to cognitive impairment (e.g., stroke, Parkinson's disease). Following autopsy, all subjects underwent neuropathological evaluation by a neuropathologist who was blinded to clinical diagnosis. Subjects received a modified (i.e., dichotomized) NIA-Reagan score based on Braak staging of neurofibrillary tangles and CERAD scoring of neuritic plaques [36,37]. In total, we identified 192 subjects with probable AD dementia due to primary AD neuropathologic change.

BPSD domains
To determine the grouping of BPSD into domains, we performed a principal component analysis with oblimin rotation on responses to BPSD questions (Supplementary Table S1) from 192 patients with dementia due to AD. We used a scree plot and Horn's Parallel Analysis to determine how many components to retain ( Supplementary Fig. S1).

Sample selection for molecular analyses
To best capture relationships between antemortem BPSD and transcriptional changes, we only considered samples from subjects with BPSD data collected within 2 years of death. Of the initial 192 subjects, we identified 100 subjects who met this criterion and had tissue available for sequencing. We created a scoring system to estimate BPSD burden for each of these 100 subjects using an additive composite of individual symptom severity and frequency within each domain. For questions with answer choices of no or yes, a score of 0 or 1 was assigned for each question, respectively. For questions with a frequency component, a low or no frequency received score of 0, a moderate frequency a score of 0.5, and a high frequency a score of 1. The individual question scores were added together to give a final score for each BPSD domain. This yielded four BPSD burden scores for each individual, one for each domain. Within each domain, we used these burden scores to identify patients who were cases (≥70 th percentile), controls (≤30 th percentile), or neither (31 st -69 th percentile), and used these groupings for subsequent transcriptome analyses. This design permitted an individual subject to be considered as a case for one domain but as a control for another, maximizing the sample sizes of cases and controls for each domain. In choosing a subset of subjects for our final sample, we made sure that each domain contained similar proportions of cases to controls and of males to females so as to minimize bias due to overrepresentation of any particular group. We ultimately chose 60 subjects for biochemical analyses (Table 1).

RNA Isolation and Sequencing
We confirmed four predominant BPSD clusters and then performed bulk RNA-seq analysis on the post-mortem ACC from a subset of 60 individuals with AD with varying burdens of BPSD. Total RNA was isolated using the QIAGEN RNeasy column-based purification kit (Germantown, MD). The quality of RNA was measured using an Agilent Bioanalyzer, which produces an RNA Integrity Number (RIN) between 1 and 10, with 10 being the highest quality samples showing the least degradation. The RINs of the 60 samples ranged between 5.3-10.0 (88% >7.0), and 1µg of high-quality RNA per sample was used for the total RNA-Seq library preparation. RNA-Seq was conducted at the Northwestern University NUSeq Core Facility. Briefly, the Illumina TruSeq Stranded Total RNA Library Preparation Kit was used to prepare sequencing libraries. The Kit procedure was performed without modifications. This procedure includes rRNA depletion, remaining RNA purification and fragmentation, cDNA synthesis, 3' end adenylation, Illumina adapter ligation, and library PCR amplification and validation. Illumina HiSeq 4000 Sequencer was used to sequence the libraries with the production of single-end 50 bp reads.

RNA-seq Differential Expression Analysis
Raw data was pre-processed with TrimGalore, including an initial quality control (QC). Read depth ranged from 50 -90M. Pseudo-alignment was performed with Kallisto [38] with k-mer 17 due to short read length (50bp). Genes were pseudoaligned to Genome Reference Consortium Human Build 38. Genes with >80% of samples with total counts <5 were removed. Principal variable component analysis (PVCA) [39] was used to identify likely important covariates by identifying factors that explained a significant proportion of variance. RIN, sex, and RNA-isolation batch were identified as contributing significantly to variation and were included in the final model, while other variables including post-mortem interval, Braak, CERAD, NIA-Reagen scores, and age at death did not. This was also analyzed with Eigen-R2 [40], which largely corroborated the PVCA conclusions. Principal component analysis (PCA) was performed and the first two principal components were visualized in scatter plot to identify likely outliers. Clustering analysis within WGCNA was also performed and were underpowered to perform differential expression analysis by sex, and therefore performed 4 independent analyses, one for each domain, using the cases and controls identified by our pre-mortem scoring system. DESeq2 [41] was used to perform differential expression analysis, and a liberal cutoff value of nominal p < 0.05 and fold-change > 0.2 was used to identify DEGs. Certain DEGs with unusually high variance and fold change were inspected for outliers, and if an isolated datapoint was >3 standard deviations from the mean, the mean was imputed for that point. Visualization of DEG overlap was facilitated with R package VennDiagram.

Functional Enrichment Analysis
Functional Enrichment Analysis was performed using gProfiler2 [42], and is described in more detail in supplementary materials.

Cellular Decomposition
Prior to cellular decomposition, WGCNA, and NetDecoder analyses, counts were converted to transcripts per million (TPM) and underwent covariate correction and variance stabilizing transformation via limma [43].
BRETIGEA [44] was then performed to estimate the relative abundance of 6 different cell types -Astrocytes, Endothelial Cells, Microglia, Neurons, Oligodendrocytes, and Oligodendrocyte Precursor Cells (OPCs)using 50 different gene markers per cell type. These results were verified with BisqueMarker [45], which agreed with did not yield significantly different trends in predicted cell type composition.

WGCNA
Weighted gene co-expression network analysis (WGCNA) was performed to identify modules of gene coexpression [46,47]. We included only the top 20% most variable genes by overall expression. Details regarding WGCNA optimization and execution are described more fully in supplementary materials. A potential hub gene was defined similar to previous guidance by WGCNA creators as having a module membership (MM) >0.8 and Gene Significance (GS) of >0.2. The notable hub genes in Fig. 3B were identified based on considerations of their high MM, high GS, overall high fold change in expression between cases and controls, and significant presence in the literature as affecting either behaviors in a particular domain or relevance to AD pathogenesis.

RNA Fluorescent Barcoding
RNA fluorescent barcoding was used to perform multiplex measurement of 37 agitation domain genes of interest (GoIs), the selection of which was informed by WGCNA and differential expression analysis. A custom CodeSet/ProbeSet (NanoString Technologies, Seattle, WA) was designed to measure GoI transcript counts from the RNA samples that remained available following use for RNAseq (N = 47). In addition to GoIs, five reference genes (IMPDH2, LAMTOR1, MTFR1L, SMIM7, TMEM50B) were selected based on low covariance between cases and controls in our RNA-seq experiment. Eight negative controls and six positive controls (NanoString Technologies) were measured as a component of QC.
Hybridization of reporter and capture probes to the RNA samples was conducted in accordance with the manufacturer's protocol (NanoString Technologies, MAN-10056-04). Briefly, 50ng of total RNA at a concentration of 10ng/µL was incubated with a Reporter CodeSet-hybridization buffer (NanoString Technologies, item no. 000136) master mix and Capture ProbeSet in a thermocycler at 65˚C for 24 hours. Incubation temperature was then reduced to 4˚C until sample processing on the following day. Hybridized samples were brought to a volume of 30µL with RNAse-free water and loaded into an nCounter SPRINT cartridge (NanoString Technologies, item no. 100078), which was run on an nCounter SPRINT profiler.
Transcript counts detected by barcode visualization in the nCounter SPRINT profiler were analyzed using nSolver Analysis Software (v4.0). Differences between group means were evaluated using Welch-Satterthwaite t-tests and a threshold of p < .05 was implemented for determination of statistical significance. All 47 measured samples were included for analysis, as binding density QC indicated sufficient RNA abundance without lane oversaturation, no fields of view were lost during imaging, and assessment of positive control linearity yielded r 2 = 1.0 for each sample.

NetDecoder
NetDecoder was performed to compare context-dependent changes in information-flow through case and control networks, as previously described [48]. A fuller description of NetDecoder is provided in supplementary materials. Briefly, we defined DEGs for each domain as source genes, and used iRefIndex v14.0 to build our interaction network. We presented the top 20 positive and top 20 negative genes in terms of flow difference or impact score for visualization across the three intermediary gene types. For visualizing changes in overall domain networks or subnetworks relating to the intermediary genes, we used Cytoscape. To simplify visualization, we filtered out edges where there was very little difference in flow, thus highlighting the results with the largest effects on the networks.

BPSD segregate into four domains
Though BPSD are heterogeneous, previous reports have indicated that common symptoms often co-occur at high rates and can be grouped into domains [23,24]. However, while groupings for BPSD are generally consistent across studies, there are slight variations that could be related to the specific cohort studied (i.e. community, nursing home, assisted living facility, etc.) or stochasticity. We performed clustering analysis of BPSD based on data from a clinical cohort where the frequency and severity of individual BPSD within 2 years from the patients' deaths was recorded. We determined grouping of BPSD into the following four domains: affective (depression and anxiety), psychosis (predominantly hallucinations), agitation (including aggression), and apathy (Fig. 1). The proportion of variance explained by each factor were comparatively similar. Delusions tended not to be explained by any single loading factor, so these were not included in the psychosis domain, although these are often grouped together clinically and in prior factor analyses [24].

Unique molecular signatures associate with each BPSD domain
After demonstrating BPSD could be split into 4 domains, we created a scoring system to estimate BPSD burden using a composite of individual symptom severity and frequency within each domain and grouped individual as cases or controls for each domain. Though BPSD are likely to result from dysfunction of multiple brain regions, the ACC has repeatedly been implicated as being involved in all four BPSD domains [49][50][51]. Therefore, we performed bulk RNA-seq on a subset of individuals from our cohort that maximized the numbers of individuals who could be considered a case or control (Supplementary Table 1 Transcriptional signatures can be used to estimate relative cell abundance of the origin tissue. Therefore, we used BRETIGEA to estimate the cell abundance of 6 major cell typesastrocytes, endothelial cells, microglia, neurons, oligodendrocytes, and oligodendrocyte precursor cells. Across all four domains, only a decrease in microglia was detected for cases compared to controls in the apathy domain, while all other cell type compositions were comparable (p < 0.05). The results were similar when a different deconvolution algorithm, BisqueMarker, was used (data not shown).
Interestingly, there were only 22 DEGs that were shared by all four BPSD domains (<3.6% of total DEGs in any domain) and all were downregulated. Functional enrichment analysis yielded pathways related to response to cytokines, integrin-mediated signaling, and collagen trimers (Fig. 1C), with TIMP1 being a notable DEG Materials for results and discussion of the other 3 BPSD domains). We found an enrichment for the extracellular matrix (ECM) including actin, collagen, glycosaminoglycans, extracellular vesicles, and cellular adhesion (Fig. 1D). Transcriptomic changes detected in agitation cases did not coincide with differences in the abundance of any individual cell type (Fig. 1E).
To confirm some of the results for our agitation domain, we quantitated transcripts using Nanostring. Though our statistical power was more limited than our initial RNAseq, we were able to confirm 17/37 genes to be significantly different between cases and controls for the agitation domain, including TIMP1, TAGLN, and FLNA (Fig. 2).
Transcriptional modules are unique and shared among BPSD domains Genes that are transcribed similarly often regulate similar biological processes, and it can be informative to group genes into co-expression modules to suggest co-regulation. In addition, these co-expression analyses can yield genes that are highly connected to the rest of the network, termed hub genes, that may be central to the transcriptional network and therefore of high interest mechanistically. We utilized WGCNA across all transcriptomes and correlated which modules were linked with the case/control condition for each BPSD domain. As expected, each BPSD domain had some modules shared across certain domains, especially for psychosis and agitation, and other modules that were uniquely significant for a single domain (Fig. 3A). In particular, an 88-gene module for growth factor and cell adhesion (greenyellow) was shared across the affective, agitation, and psychosis domains; a 28-gene module for ECM and actin cytoskeleton (darkturquoise) was shared across apathy, agitation, and psychosis; and a 98-gene module for factor activity (magenta) was shared across agitation and psychosis; no module was associated with all four domains, again suggesting the separability of these traits on a transcriptional level. All shared modules suggested a reduction in transcription in cases, consistent with differential expression trends. Focusing on the agitation domain, three modules were uniquely correlated: modules for nucleosome assembly (lightyellow; 35 genes), ATPase and synaptic signaling (purple; 916 genes), and post-synaptic signaling and response to monoamines (darkgrey; 63-gene). The postsynaptic and monoamine module was enriched for serotoninergic signaling, dopaminergic signaling, and GABAergic signaling. Interestingly, the synaptic signaling modules (purple and darkgrey) were amongst the only significantly correlated modules that demonstrated increased transcription for cases.
Hub genes in WGCNA that are connected with a given domain are defined by the high connectivity within their module as well as significant influence on transcription with the BPSD domain (Fig. 3B). The ECM module, associated strongly with agitation and psychosis and more weakly with apathy, yielded 10 potential hub genes, including TAGLN (MFC,agitation = 0.362, p = 1.3 x 10 -5 ) and FLNA (MFC, agitation = 0.646, p = 5.9 x 10 -6 ). TIMP1 was also in this domain but had module membership value slightly below the cut-off for potential hub genes.

Potential drivers of transcriptional network information flow with agitation domain
BPSD likely arise from complex changes in multiple biological networks, and the dynamic interactions of genes and proteins within a network may be key to understanding the difference between the BPSD cases and controls in AD. While enrichment analyses suggest important pathways based on overrepresentation of genes compared to chance, they cannot inform how biological information may change in a context-dependent mannersuch as a disease vs non-disease state. Therefore, we used NetDecoder to compare the 'information-flow' through each BPSD domain's cases and controls' networks, which utilizes a process-guided flow algorithm to identify the weights of information flow from source genes, DEGs, to target genes, transcriptional regulators [48]. For context, when this algorithm was applied to an older breast cancer dataset, they were able to identify three laterconfirmed, prognostic markers as important drivers of information flow in the cancer network that were not originally implicated in the original study that lacked a context dependent approach. Important genes in NetDecoder are labeled as high impact genes, network routers, or key targets, and by function of the algorithm, all are genes that were not identified as DEGs but likely affect overall information-flow through the biological system via regulation of downstream transcription.
Focusing on the agitation domain, NetDecoder revealed divergent information flow between cases and controls, and the top 40 network routers, key targets, and high impact genes are shown (Fig. 4 A-E; Supplementary Fig.   S3). Though a number of high impact genes are of interest, two are particularly notable (Fig. 4 E, Arrows   indicated). The ER-beta, encoded by ESR1, was identified as the top key target and high impact gene mediating positive information flow (Fig. 4F). While little remains known about ESR1 function in agitation/aggression in the frontal cortex, ESR1 + cells within the hypothalamus and other limbic regions have deep evidence linking them to control of aggressive behaviors [52][53][54][55][56][57]. The other is Parkin, the ubiquitin ligase encoded by PARK2 (Fig. 4G) and genetically associated with familial Parkinson's disease. PARK2 has significant cross-talk with tau, regulates mitophagy in AD, and interestingly has been linked to impulsive behaviors in Parkinson's disease [58][59][60].

DISCUSSION
The data presented here represents the first exploration of affective, apathy, and agitation symptoms on the transcriptome-wide level and establishes unique patterns of mRNA expression in one of the most consistently implicated brain regions to BPSD, the ACC [49][50][51]. In addition, we add to the growing knowledge base about transcriptional changes in AD with psychosis [28] (discussed in greater detail along with the affective and apathy domain in supplementary materials). We confirmed that commonly co-occurring BPSD symptoms cluster into domains in our cohort and that these domains are typified by unique transcriptional signatures in the ACC, even when individual samples are used in an overlapping design. Using co-expression analysis, we observed individual BPSD domains being associated with shared and unique transcriptional modules with potential hub genes that may serve as targets for future discovery of discrete molecular mechanisms and novel pharmacology. Finally, we identified key drivers of information flow through biological networks associated with BPSD domains, highlighted by ESR1 and PARK2 being potential mediators of the agitation domain in AD.
Though BPSD may be thought of as manifestations of late life primary psychiatric disorders (PPD), there is already some evidence suggesting that BPSD mechanisms diverge molecularly. There have been a handful of genetic studies suggesting overlapping polygenic risk for PPD and neurodegenerative disease [61,62], and while good epidemiological evidence suggests PPD are risk factors for resultant neurodegenerative dementias [63,64], the few studies comparing PPD and BPSD on a genetic level have yielded surprising results. For instance, while psychosis -focusing on the 'positive' symptoms of hallucinations and delusionsis a prominent symptom in schizophrenia, bipolar disorder, and AD with psychosis, a recent GWAS found negative correlations between schizophrenia and AD with psychosis and bipolar disorder and AD with psychosis, suggesting not just a lack of association between these PPD and BPSD, but a reduced risk of psychosis in AD with increased polygenic risk for schizophrenia or bipolar disorder [65,66]. This is a stark contrast to the extensive genetic overlap between schizophrenia and bipolar disorder risk [67], suggesting that BPSD may mechanistically distinct from PPD.
Treatment of BPSD has also yielded surprising differences from PPD. For instance, the HTA-SADD trial found no benefit of two common serotonergic antidepressants for depression in AD [68], and a Cochrane Database meta-analysis support this lack of efficacy [69]. Similarly, while selective serotonin reuptake inhibitors (SSRIs) are not used to treat hallucinations and delusions in PPD for psychotic disorders, and may even induce psychosis in bipolar disorder through increasing the risk of mania [70], a common SSRI citalopram seems to have some benefit for reducing these psychotic symptoms in AD, though this was discovered on secondary analysis and requires follow-up[71]. Similar to genetic differences between BPSD and PPD, it seems likely that pharmaceutical approaches need to differ substantially in treating the two disorders, which necessitates the further investigation of BPSD as its own entity distinct from PPD.
Rigorous characterizations of BPSD antemortem with complementary post-mortem tissues for molecular analyses are remarkably scarce, so we sought to optimize our analytical power for the samples we were able to obtain. In addition, given the high prevalence of BPSD in AD ( > 95%), finding enough controls without any BPSD would be very challenging. This required us to adopt an experimental design where each individual's sample could be considered to be a case (25% highest score of the domain) or control (25% lowest score of the domain) for a specific BPSD domain, meaning that some transcriptome could be potentially analyzed multiple times depending on the comparison. Therefore, it was encouraging to see such a large number of DEGs that are unique to each domain despite the overlapping design, which may suggest that each BPSD domain has distinct biological etiologies despite the common neuropathological driversin this case, AD. These findings may be comparable to other investigations showing distinct pathways associated with cognitive decline in AD at the transcriptomic, proteomic, and methylation levels in brain tissue without overlapping dependence on being correlated with AD pathology [72][73][74][75][76][77][78][79], highlighting the heterogeneity of downstream molecular processes from what is putatively considered the upstream etiological agents, namely A and tau. Though this design has the advantage of increasing statistical power in exploring a wide array of symptoms, it presumes separability that precludes identification of individuals with overlapping psychiatric domains that may have unique molecular mechanisms distinct from if these symptoms presented separately. Future investigations can help clarify this.
Given this overlapping design, it was surprising that so few genes were shared amongst all the BPSD domains.
It is notable that response to cytokines was implicated as a functionally enriched pathway, as neuroinflammation is often considered an integral driver of neurodegeneration and subsequent synaptic dysfunction [80]. TIMP1 was among the 22 shared DEGs, all of which were downregulated. A major inhibitor to a number of metalloproteases such as MMP3 and MMP9, TIMP1 has been implicated in multiple forms of neurodegeneration and neuroinflammation [81,82], with the hypothesis that early upregulation of TIMP1 maintains balance in neurodegenerative states while late downregulation may suggest inability to achieve homeostasis [83]. It is therefore speculative but possible that reduced expression of TIMP1 leads to loss of homeostasis, which could lead to heightened stochasticity in downstream processes and divergent BPSD. While it is also interesting that TIMP1 has been suggested as a biomarker in biofluids in both Parkinson's disease and AD [84][85][86], further exploration of TIMP1 could be especially fruitful in understanding the early drivers of BPSD.
Though our analyses discovered multiple interesting targets that were unique to each BPSD domain (see supplementary materials for in depth results and discussion), we focused on the results of the agitation domain.
The enrichment analysis of the DEGs for this domain were highly suggestive of changes in the ECM or matrisome, and a shared transcriptional module associated with agitation, apathy, and psychosis similarly was enriched for the ECM. A recent and extensive proteome-wide study in AD found that a module of co-expressed proteins enriched for the matrisome was highly correlated with global AD pathology and ApoE status but shockingly independent of cognition [87]. The possibility exists that changes in the ECM are less correlated with cognition but are better tied to BPSD, especially agitation. Interestingly, a common functional SNP in MMP9, which is inhibited by TIMP1, was associated with inhibition of aggression and irritability in one study [88].
Currently, much more evidence would be needed to link changes in the ECM with agitation in AD, but this finding further highlights the necessity of considering BPSD in studies of neurodegeneration in addition to cognition, as there may be diverging molecular mechanism related to both sets of symptoms.
The a module for post-synaptic signaling and monoamines was upregulated in agitation cases versus controls and was enriched genes that are common treatment targets for agitation, such as those related to dopamine (antipsychotics), serotonin (antipsychotics and antidepressants), and GABA-A receptors (benzodiazepines) [22].
As monoaminergic treatments are so important in PPD but have limited to no success for the affective domains [68,69], it was interesting to note that none of the serotonergic, dopaminergic, adrenergic, or muscarinic receptors were DEGs for the affective domain. In contrast, DRD1 and 5HTR2C were DEGs for the psychosis and agitation domains while the psychosis domain also demonstrated DRD2 and DRD4 as DEGs.
While far from conclusive, this again dovetails with the treatment failures for SSRIs for affective behavior in BPSD and further supports the hypothesis that PPD have fewer mechanistic similarities to BPSD.
Additionally, the finding that DRD1, a GPCR coupled to Gs/, and downstream effectors PDE1B and CAMK4 are associated with this module and the agitation domain is in line with some genetic reports that DRD1 SNPs are associated with greater impulsivity and aggression [89,90], including in AD [91,92] and Parkinson's disease [93]. The role of striatal role of DRD1 in aggression has been demonstrated before [94], but how DRD1 and some of its downstream signaling molecules affect agitation/aggression in the ACC is less clear. Similar to DRD1, increased GABA-A signaling in the prefrontal cortex, including the ACC, has been linked with increased aggressive and impulsive behaviors [95,96] and may interact with CAMK4 expression in certain situations [97]. It was notable that those with significant agitation domain behaviors trended towards having more neurons in the ACC than controls, though this analysis precluded investigation of the neuronal subtypes. It is possible that relative increases in GABAergic neurons in cases or decreases in glutamatergic neurons in controls could lead to these differences in agitation domain behavior.
Transcriptome datasets are inherently noisy, and predicting how biological networks regulating transcription, translation, and protein-protein interactions might change based on differential gene expression can uncover hidden drivers of disease. Using an analytical framework to analyze transcription in this context-dependent way, we uncovered two high interest gene targets related to the agitation domain, ESR1 and PARK2. For neurons in the ventromedial hypothalamus [52], posterior and medial amygdala [56], and bed nucleus of the stria terminalis [57], ESR1 + expression differentiates these cells as ones that regulate aggression from the behavioral function of ESR1cells. Despite the importance of ESR1 as a marker of aggression-regulating neurons, the actual function of ESR1 in the cell and resultant aggressive behavior is less clear, though knockout of ESR1 in mice leads to reduced aggression [98] and ESR1 polymorphisms have been linked with aggression in humans [99] and songbirds [100]. Understanding ESR1's role in the ACC in terms of agitation, aggression, and impulsivity may be a particularly fruitful avenue for mechanistic understanding and future drug development.
The main novelty of our study is that it is the first, to our knowledge, to study the affective, apathy, and agitation domains on the transcriptome-wide level, which is comparable to the many studies of cognition in AD [73,75,77], dating back at least as early as 2008 [72], and the recent study of AD with psychosis [28].
Despite this, there are important limitations worth noting. First, our bulk tissue approach precludes deeper exploration of the role of different cell types in BPSD. Similarly, we present only one brain region's worth of data, and it is likely that interactions with other regions are necessary to fully understand each BPSD domain's pathogenesis. We were also unable to correlate our findings with medication data before the patient's death, which will be an important covariate to include in future studies. Additionally, transcriptomic differences do not always translate into protein differences [87], and so future proteomic studies would help solidify the significance of our results. Even with these limitations in mind, we hope our work will lead to future molecular investigations into BPSD so that advanced therapeutics can be designed and translated to the clinic for many of our society's most vulnerable and affected patients.  Fig. 2 Validation of key DEGs for the agitation domain. The differential expression of A) Six DEGs that are also potential hub genes in the ECM (Darkturquoise) module and B) Five DEGs that are potential hub genes in the post-synapse (Darkgrey) module were confirmed with RNA fluorescent barcoding after expression was normalized to five housekeeping genes. C) Non-WGCNA genes and modules with a single confirmed potential hub gene (i.e., ADIRF from the nucleosome (Lightyellow) module and CD33 from the transcription factor (Magenta) module) are represented. Blue circle = male, red circle = female; **** = p < 0.0001, *** = p < 0.001, ** = p <0 .01, * = p < 0.05.